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DFTS  ON  IRREGULAR  GRIDS:  THE  ANTERPOLATED  DFT 

VAN  EM DEN   HENSON  * 

Abstract.  In  many  instances  the  discrete  Fourier  transform  (DFT)  is  desired  for  a  data  set 
that  occurs  on  an  irregular  grid.  Commonly  the  data  are  interpolated  to  a  regular  grid,  and  a  fast 
Fourier  transform  f  F FT)  is  then  applied.  A  drawback  to  this  approach  is  that  typically  the  data  have 
unknown  smoothness  properties,  so  that  the  error  in  the  interpolation  is  unknown. 

An  alternative  method  is  presented,  based  upon  multilevel  integration  techniques  introduced  by 
A.  Brandt.  In  this  approach,  the  kernel,  e-"*',  is  interpolated  to  the  irregular  grid,  rather  than 
interpolating  the  data  to  the  regular  grid.  This  may  be  accomplished  by  pre-multiplying  the  data  by 
the  adjoint  of  the  interpolation  matrix  (a  process  dubbed  anterpolation),  producing  a  new  regular- 
grid  function,  and  then  applying  a  standard  FFT  to  the  new  function.  Since  the  kernel  is  C'°°  the 
operation  may  be  carried  out  to  any  preselected  accuracy. 

A  simple  optimization  problem  can  be  solved  to  select  the  problem  parameters  in  an  efficient  way. 
If  the  requirements  of  accuracy  are  not  strict,  or  if  a  small  bandwidth  is  of  interest,  the  method  can 
be  used  in  place  of  an  FFT  even  when  the  data  are  regularly  spaced. 

1.   The  formal  DFT  and  the  ADFT.   The  DFT  is  defined  as  an  operation 

that  maps  a  Iength-iV  complex- valued  sequence  { ru.  X] -r.v-i }  to  another  length- 

.V  complex-valued  sequence  {}o.}\ '".v-i}  by  the  rule 

,v-i 

(i)  ik  =  Y,  V~,J~A/'V-    for  fc  =  o,i v-i. 

As  defined  in  (1).  the  DFT  is  performed  on  data  that  are  presumed  to  be  given 
on  a  regular  grid,  with  constant  spacing  between  the  data  points.    Furthermore,  the 

transform  values  {io-^i J'.v-i}  are  also  presumed  to  lie  on  a  regular  grid  in  the 

frequency  domain.  In  many  applications,  however,  the  data  for  a  problem  are  not 
spaced  regularly.  It  is  of  some  interest,  then,  to  determine  how  a  discrete  Fourier 
transform  may  be  computed  for  such  a  data  set.  To  perform  this  computation,  we 
develop  and  implement  in  one  dimension  an  algorithm  based  on  multilevel  integration 
techniques  outlined  by  Achi  Brandt  ([2].  [1]).  The  method  presented  here  can  also  be 
developed  for  higher-dimensional  problems.  One  application  of  this  technique  [G]  is  in 
the  reconstruction  of  images  from  projections  (inverting  the  Radon  transform). 

To  begin,  it  is  necessary  to  decide  what  is  meant  by  a  Discrete  Fourier  Transform 
for  irregularly  spaced  data.  Therefore,  the  concept  of  a  formal  DFT  is  introduced, 
which  is  defined  as  follows: 

Consider  any  set  of  N  ordered  points  in  the  interval  [O.A').  satisfying 

0  <  Jo  <  Ji  •  •  •  <  J.v-i  <  -V 

and  suppose  a  vector- valued  function  (grid  function)  u(Xj)  is  specified.  The  formal 
DFT  is  defined  as  the  M  =  M\  +  M2  +  1  quantities 

(2)  u(w,)  =  ]T  u{xj)e-^x>  ,       w,  -  — ,    -Af,  <  /  <  M2  . 

j=o  A 

where  /  is  an  integer,  and  M\  and  M2  are  positive  integers  specifying  the  range  of 
frequencies  of  interest.  The  formal  DFT  may  be  thought  of  as  an  approximation  to  a 
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selection  (  —  Mi  <  /  <  A/2)  of  the  Fourier  coefficients  of  u(x).  In  this  view,  u(x)  should 
be  regarded  as  an  A'-periodic  function  known  only  at  the  grid  points  x;. 

It  is  desired  that  this  sum  be  calculated  to  a  prescribed  accuracy,  say  f||u||i,  where 
||u||i  is  the  discrete  L\  norm  ||u||i  =  A-1  J2j^  lu(zj)l-  Note  that  any  grid  spacing 
is  allowed  for  the  x}  (in  particular,  the  spacing  needn't  be  constant),  that  there  is  no 
relationship  between  the  integers  Afi,A/2,  and  N,  and  that  there  is  no  requirement 
that  these  integers  have  any  special  value,  (such  as  being  powers  of  2).  Calculating  the 
sum  in  (2)  directly  would  have  a  computational  cost  of  O(MN)  operations.  Instead 
of  forming  this  sum  directly,  though,  an  approximation  to  it  will  be  computed,  using 
an  F FT  to  accelerate  the  computation. 

The  procedure  begins  with  the  definition  of  an  auxiliary  grid,  fi;y  ,  covering  the 
range  of  values  [0,A').  Let  A.  be  an  integer,  whose  value  will  be  determined  shortly, 
and  let  the  grid  spacing  h  be  defined  by  h  =  X/N+.  Then  the  auxiliary  grid  consists 
of  the  points  yT  -  (r  —  l)/i,  for  r  =  1,2, ...  A',. 

Suppose  that  the  value  of  some  function  g(yT)  on  the  regular  grid  ftjy  is  to  be 
interpolated  to  the  gridpoints  Xj  by  Lagrangian  interpolation.  We  will  identify  the 
interpolation  by  specifying  the  degree  of  the  polynomial  to  be  used.  Thus,  p-degree 
Lagrangian  interpolation  is  computed  using  a  polynomial  of  degree  p  or  less.  For  each 
Xj,  the  p  +  1  nearest  neighbors  on  the  grid  fl1^-  must  be  located.  Let  these  points  be 
designated  yK(J<o).  y^(jj)-,  ■  ■■  ,yK(j,p)-  These  points  should  be  chosen  modulo  X,  so  that 
a  point  near  the  limits  of  the  interval  [0,  A')  may  have  neighbors  near  the  other  end 
of  the  interval.  (This  is  justified  because,  as  will  be  seen  shortly,  the  function  to  be 
interpolated  for  the  formal  DFT  is  A'-periodic.)  For  each  Xj,  the  p  +  1  Lagrangian 
interpolation  weights  are  computed  by 

(3)  wB(xi)=       [  ; :    • 

m=o  {**&■*)  ~  y*u,m)> 

and  the  interpolation  of  the  function  g  to  the  gridpoints  Xj  is  given 

v 
g(xj)  %  ^wn(xj)y(yK(iin))  . 

n=0 

Letting  g  be  the  vector  of  function  values  g{>jk)  and  g  be  the  vector  of  interpolated 
values  g(xj),  the  interpolation  may  be  written  in  matrix  form 

9  =  [iy]g 

where  1*  is  the  N  x  A,  interpolation  matrix  mapping  a  function  on  H^  to  the 
gridpoints  {xj}.  The  entries  of  this  matrix  are 

lJT]       _  J  wn(xj)     if  n(j,n)  =  m 
1  »]jm  ~  \  0  else   . 

We  are  now  ready  to  compute  an  approximation  to  equation  (2).  The  strategy 
will  be  to  interpolate,  for  each  u/,  values  of  the  kernel  e~XuJ[I)  from  the  auxiliary  grid 
ftyv^.  That  is,  equation  (2)  is  approximated  by 

N-\ 

(4)  u(w/)»  Yl  u(xj)e("l,Xj) 

j=o 

2 


where 

p 
e(uhXj)=  Trwn(Xj)e-iwtv<>rt  . 

n  =  0 

Let  the  column  vector  of  exponential  values  e~,UJ,Vk  be  designated  c/,  and  the  vector 
of  interpolated  kernel  values  e(u>i,Xj)  be  denoted  t\  Then  this  interpolation  can  be 
written 

e,  =  [lxy]e,   . 

Notice,  however,  that  c/  is  to  be  used  as  the  ltfl  row  of  the  matrix  giving  the  kernel 
of  the  summation  in  (1).  To  compute  t\  as  a  row  vector,  the  adjoint  of  the  matrix 
equation  equation  is  needed,  namely 

(ei)T  =  (2J?/)T  =  («i)T[2J]T  • 

Lot  us  define  the  Af-vector  u  =  u(w/),  the  N-vector  u  =  u{t ,).  and  t ho  matrices  giving 
the  kernels  of  equations  (2)  and  (4)  as  W  and  U'.  respectively.  Let  the  M  x  .V.  matrix 
whose  Z''1  row  consists  of  e_lLj'y"  for  the  A',  points  on  fly  be  designated  IF  (it  is  useful 
to  observe  that  this  matrix  consists  of  M  consecutive  rows  of  the  standard  1) FT  kernel 
for  a  uniform  grid  with  A',  points).   Then  the  formal  DFT  is  approximated 

u     —     \\  u 
as     Wu 

=  mmTu . 

This  notation  can  be  simplified  slightly  by  denoting  the  vector  created  by  multiplying 
f7  by  [IlY ' .  as  u.  Since  the  matrix  [I^]T  is  the  adjoint  of  the  Lagrangian  interpolation 
matrix,  the  process  of  computing  u  —  [2*]  u  has  been  dubbed  (interpolation.  Then 
the  approximation  to  the  formal  DFT  is 

( 5 )  u  %  W  u   . 

which  we  call  the  Anterpolated  Discrete  Fourier  Transform  (A DFT). 

The  A  DFT.  as  a  matrix  multiplication,  requires  0(MXm)  operations.  In  general, 
A'»  will  exceed  Ar,  so  as  a  matrix- vector  multiplication,  the  AD  FT  has  no  advantage 
over  (2).  If,  however,  Ar,  is  selected  appropriately,  the  approximation  can  be  computed 
quite  rapidly.  Let  M,  =  max  {Mi,  A/2}.  Then  if  N»  is  selected  such  that  A'.  >  2A/», 
and  at  the  same  time  JV„  is  a  number  for  which  an  FFT  module  exists,  then  the  fast 
Fourier  transform  can  be  applied  to  compute  the  DFT  summation 

FFT  {u}  =  ir 'fiTe-*"W    for  /=-— +1,-— +  2,...—   . 
A.   rt-0  2  2  2 

Recalling  that  h  —  A'/ A*,,  it  may  be  seen  that  the  DFT  summation  therefore  yields 
( l/A'.)u((27r/)/A' ).  Multiplying  by  AT,  thus  yields  a  set  of  values  that  includes,  as  a 
subset,  all  the  desired  values  of  u(u>i). 

Computing  the  ADFT,  then,  consists  of  two  phases: 

1  -   u  is  computed  from  u  by  anterpolation:  u  =  [2y]Tu. 

2.  u  is  computed  from  u  by  a  Fast  Fourier  Transform. 
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2.  Operation  count  for  the  ADFT.  The  cost  of  computing  the  ADFT  con- 
sists of  the  cost  of  computing  the  interpolation  weights,  the  cost  of  computing  the 
vector  u  —  [I* }Tu.  and  the  cost  of  the  FFT  on  N*  points. 

Computing  the  (p  +  1)N  interpolation  weights,  wn(xj),  by  the  formula  in  (3)  is 
the  cost  of  the  computing  the  numerator,  since  the  regular  spacing  on  fi^  means 

that  the  denominators  of  wn(xj)  are  independent  of  j.  To  compute  the  numerators, 

p 
the  product      T  (xj  —  xK/jiTn\)  is  computed  for  each  Zj,  requiring  2p  +  1  operations. 

m-0 

Then  the  ntk  interpolation  weight  can  be  obtained  by  dividing  by  the  product  of 

p 

(xm  -  xK(jtTn))  with  the  precomputed  denominator      J  (xK(nj\  —  xK^j^m\),  requiring  2 

m  =  0 
m  ^n 

operations  for  each  of  the  p  +  1  weights  associated  with  the  point  Xj.  The  calculation 
of  the  weights  thus  requires  0(N{p+  1))  operations.  It  is  important,  however,  to  note 
that  the  calculation  of  the  weights  is  dependent  only  on  the  relationships  between  the 
gridpoints  {y/,.}  and  {xj},  and  is  independent  of  the  data  set,  u(xj).  This  means  that 
if  a  known  set  of  gridpoints  {x_,}  and  a  standard  auxiliary  grid  Sly  are  to  be  used 
repeatedly,  the  interpolation  weights  wn(zj)  may  be  precomputed  and  stored,  and 
needn't  be  included  in  the  cost  of  the  algorithm.  This  will  be  assumed  to  be  the  case. 
The  matrix  [J5]  is  A",  x  .V  and  the  data  vector  il  is  .V  x  1.  so  the  computation  of 
u  =  [Iy}Tu  would  be  0(\ N.)  if  performed  as  a  matrix-vector  multiplication.  There 
is,  however,  a  much  more  efficient  method.  The  index  table  k(j.ti)  can  be  stored 
along  with  the  interpolation  weights.  For  each  Xj  and  for  each  /;.  the  value  of  k(j,  n) 
is  the  index  of  the  nth  interpolation  neighbor  that  is  used  to  interpolate  from  Slv< 
to  the  gridpoint  Xj.  The  periodic  nature  of  the  kernel  being  interpolated  means  that 
the  interpolation  is  always  to  a  gridpoint  x^  in  the  center  of  the  set  of  p  interpolation 
points  (as  is  well  known.  [5],  the  Lagrangian  interpolation  is  better  behaved  when  this 
is  the  case).  If  p  is  odd,  then  x,  alwavs  lies  between  y  ,    p-i ,  and  u  ,    „+\ .,  while  if  p 

is  even,  then  >JK(j_t)  is  the  closest  gridpoint  on  f]v<  to  Xj. 

Computing  the  vector  u  is  then  very  easy,  and  may  be  done  in  2N(p+  1  j  operations, 
according  to  the  algorithm: 

1.  Initialize  u(yT)  -  0  for  all  yT  G  fthv>.  (1  <  r  <  .V.) 

2.  For  j  =  0,1,...,N-  1 

For  n  =  0, 1, .. .  ,p 

Set  u(yK(j>n))    —     u(yK(j,n))  +  ™n{Xj)u(xj). 

Having  computed  the  values  of  it,  consider  now  the  cost  of  the  FFT  portion  of 
the  ADFT.  This  is  simply  the  cost  of  an  A'.-point  FFT.  In  the  next  section  criteria 
for  choosing  Nm  will  be  determined.  For  now  the  only  requirement  is  that  an  FFT  can 
be  computed  on  a  vector  of  length  A7..  As  such,  N.  must  have  factors  for  which  FFT 
modules  are  available.  For  the  purpose  of  an  operation  count,  however,  it  is  easiest 
to  assume  that  A,  is  a  power  of  2.  Indeed,  we  shall  see  that  we  have  great  flexibility 
in  our  selection  of  Ar»,  and  since  powers  of  2  or  4  produce  the  most  efficient  FFTs, 
this  is  a  good  assumption.  In  this  case,  the  cost  of  the  FFT  portion  of  the  ADFT  is 
0(A.  log2A,). 


The  costs  of  the  ADFT  can  now  be  computed.  In  terms  of  data  storage  it  requires 
four  arrays.  One  is  the  A-point  vector  containing  the  input  data.  u.  In  addition,  an 
A.-point  complex  vector  is  required  for  the  input  and  output  of  the  F FT .  Assuming 
that  the  weights  are  precomputed  and  stored,  two  auxiliary  arrays  are  necessary,  an 
A  X  (p  +  1)  real  (or  double  precision)  array  holding  the  interpolation  weights  wn(ij), 
and  the  A  X  (p+  1)  integer  array  of  indices,  k(j,ti). 

If  the  operations  of  multiplication  and  addition  are  counted  equally,  and  if  the 
weights  and  indices  are  pre-stored,  the  operation  count  is  C\S.  log2  A",  complex  oper- 
ations for  the  F FT  portion  of  the  algorithm,  where  C\  depends  on  the  choice  of  FFT 
algorithm.  The  computation  of  u  entails  2\(p  +  1 )  operations  that  are  real  or  complex 
according  to  whether  u  is  real  or  complex.  Counting  both  phases  of  the  algorithm, 
the  operation  count  of  the  AD  FT  is 

CiN.  log2A.  +  2N(p+l)  . 

This  should  be  compared  with  the  operation  count  of  the  formal  DFT,  which  is 
0(M A).  The  computation  of  the  ADFT  is  more  efficient  provided  M  is  larger  than 
2(p+  1)  +  ( A'„  log2  A.  )/.Y.  a  condition  that  will  generally  occur  in  practice. 

3.  Error  analysis  for  the  Anterpolated  DFT.  One  of  the  attractive  fea- 
tures of  the  ADFT  is  that  the  interpolation  is  performed  on  the  kernel,  which  has 
known  smoothness  properties,  rather  than  the  data  set.  which  generally  has  unknown 
smoothness  properties.  Since  interpolation  error  depends  on  the  smoothness  of  the 
interpolated  function,  the  error  committed  by  using  the  ADFT  is  relatively  easy  to 
analyse. 

Consider  the  error  in  /j-degree  Lagrangian  polynomial  interpolation,  when  the 
interpolation  is  from  a  set  of  p  +  1  gridpoints  that  are  equally  spaced.     Let   these 

gridpoints  be  designated  so.£i s;>-    A  function  /(x),  whose  values  are  known  at 

these  gridpoints.  is  to  be  interpolated  to  the  point  x  £  [sO-s>]-  Let  x  =  £o  +  th,  where 
h  is  the  gridspacing,  and  /  £  [O.p).  The  approximation  to  /(£o  +  th)  's  the  value  of 
the  Lagrangian  interpolation  polynomial  Pp(£o  +  th), 

PP(z)  =  f>i(*)/(fc)     where     w,-(x)  =   f[  \*  "  ^  \  ■ 

m  ^  I 

Defining  7rp(r)  =  t(t  -  1  ){t  -  2) .  . .  (/  -p),  and  (  £  [^o<^p],  the  error  in  the  interpolation 
is  bounded  by 

(6)  |/(6  +  th)  -  PP($o  +  th)\  <  |Tp(/;ll^+1l'9p+l 

(P+  1): 

where  Qp+\  =  max^^  ]|/^P+1)(C)I-  See  [8],  pages  264-270,  for  a  derivation  of  this 
error  term. 

It  is  useful  to  bound  this  error  more  precisely.  To  do  this,  we  examine  the  behavior 
of  the  factorial  polynomial  np{t).  This  polynomial  has  been  well-studied,  and  many 
results  can  be  found  in  various  numerical  analysis  texts,  ([5],  [S],  [9],  [11]).  These 
results,  however,  are  developed  for  the  case  that  x  can  be  anywhere  in  [£o,fp].  In  the 
present  case  the  interpolation  is  always  to  the  center  subinterval.  Thus  for  p  odd, 
t  £  [^,^1,  while  for  p  even,  either  t  £  [§,§+  1]  or  t  €  [f  -  l,f]. 
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To  shorten  the  discussion,  assume  that  p  is  odd.  This  is  the  most  common  case, 
where  p  =  1  gives  linear  interpolation  and  p  =  3  give  cubic  interpolation.  Similar 
results  can  be  obtained  for  p  even.  Consider  the  following  lemma,  the  proof  of  which 
may  be  found  in  [6]. 

Lemma  1.   If  p  is  a  positive  odd  integer,  then 


te[iHi,£±i]  (p  +  1)! 


max 


< 


1 


2p+i 


This  result  can  be  used  to  find  an  error  bound  for  the  ADFT.   In  this  case  the 
functions  being  interpolated  are 

.V  Y  V 

e~l^x   for  /=  -—  +  1,-—  +  2,...—     . 
2  2  2 

From  this  set  of  functions,  the  only  ones  whose  values  are  of  interest  are  those  for 
/  between  —Mi  and  A^.  Recalling  the  definition  M.  —  max  {M\ .  A/2},  the  largest 
absolute  value  among  the  frequencies  of  interest  is  u,',\/  =  (2irMm)/X .  Therefore 


max 

■r6[yK(j.o)-yKO,pj] 


rp+l  I 


9 


l^A/ 


iP+1 


Insprting  this  and  the  bound  from  Lemma  1  into  equation  (6)  gives 

,  /iw\f  \p+1 

(7)  |e-<^-eU.*,)|     < 

which  is  then  used  to  obtain 

\u-Wu\     =     IEjL'o1  u(xj)e 


v-i  ..  _IW|Ij 


Ejlo1  "Uj  m~/. -01 


=     (/i^A//2)f+1.V||u||1   . 

Finally,  substituting  /?  =  X/iV»  and  wm  =  (2ttM*)/X  establishes  the  desired  error 
bound.  The  error  in  the  ,4 DFT  approximation  to  the  formal  DFT  is  bounded,  for 
-Mx  <  I  <  M2.  by 


(8) 


l^r)P  ljv|H|1, 


where  w/  =  2tt//A\  and  the  ADFT  (5)  is  computed  using  an  FFT  of  length  Ar, 
Since  the  bound  holds  for  all  desired  values  of  W/,  it  then  follows  that 


(9) 


l«-[^]|L<(^)P+1^IHIi  , 


where  II  •  II      is  defined  as  the  maximum  absolute  value  in  the  vector.  It  is  also  worth 

II         II  oc 

noting  that  an  error  bound  for  any  desired  frequency  can  be  obtained  by  replacing 
cj,v/  with  ui  in  the  derivation,  leading  to 

10] 


uU;)-[\Vu}(^)\<  (j^j         N\\u\U 


This  is  especially  useful  information  for  those  occasions  when  only  the  low  frequency 
components  are  of  interest,  or  when  the  accuracy  required  of  the  approximation  is 
greater  for  the  low  frequencies  than  for  the  high  frequencies. 

4.  Selection  of  p  and  .V..  The  error  bound  just  derived  is  useful  in  that  it 
provides  a  way  to  select  the  operational  parameters  X.  and  p.  Recall  that  the  goal 
is  to  calculate  an  approximation  to  u  to  some  prescribed  accuracy,  |u  -  Wu\  <  (\\u\\. 
In  practice  we  will  want  to  make  the  error  small,  so  it  will  be  assumed  that  e  <C  1. 
Comparison  with  (8)  gives  the  requirement 

A/.t\''+i        e 
< 


.  .V.    }  -  X 

which  may  be  written  as 


(11)  .V.  >  .\L- 


.V\5±r 


i 


For  a  given  formal  DFT,  the  values  of  A/„.  .V.  and  <  are  considered  to  be  part  of 
the  problem  specification.  To  ensure  that  the  specified  accuracy  is  obtained,  it  is  only 
necessary  to  select  integers  jV»  and  p  so  that  (11)  is  satisfied.  Naturally,  there  may  be 
many  combinations  of  parameter  values  that  achieve  this  goal.  The  parameters  should 
therefore  be  selected  to  fullfill  some  other  desirable  property  as  well.  Specifically,  they 
should  be  selected  also  to  minimize  the  computational  effort  of  the  algorithm. 

To  see  how  this  may  be  accomplished,  recall  that  the  work  involved  in  comput- 
ing the  ADFT  with  .V,  points  on  the  auxiliary  grid  ftjy  and  p-degree  Lagrangian 
interpolation  is 

0(Nm\og2Nm)  +  0{N{p+  1))  . 

The  value  of  the  constant  on  the  0(N(p+  1))  term  depends  on  whether  the  weights 
and  indices  are  pre-stored,  or  calculated  "on  the  fly".  For  the  analysis  that  follows, 
we  assume  the  weights  and  indices  are  pre-stored.  in  which  case  the  constant  is  2. 
The  constant  on  the  first  term  depends  on  several  factors.  FFTs  generally  have 
a  complexity  of  {X /q)\og{X/q)  for  some  number  q  >  1.  If  the  data  have  certain 
symmetries,  then  a  specialized  FFT  may  be  used  for  faster  computation  ([3],  [7], [10]). 
The  variety  of  available  FFT  algorithms  pursuades  us  to  leave  the  constant  on  the 
first  term  as  an  unspecified  parameter,  C\. 

The  total  work  in  computing  the  ADFT  can  therefore  be  written  as  a  function  of 
the  two  parameters  X,  and  p.  For  a  fixed  problem  size  (X  and  A/,),  and  a  prescribed 
error  tolerance  c,  the  work  in  computing  the  ADFT  to  the  required  accuracy  is 

(12)  \V(X.,p)  =  C1X.\og2X.  +  2X(p+\)   , 
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and  we  seek  an  optimal  parameters  minimizing  W(Nm,p)  over  all  combinations  (N+,p) 
satisfying  (11),  if  such  a  choice  exists. 

Limiting  cases  may  be  determined  by  examining  nearest  neighbor  interpolation 
(p  =  0),  as  well  as  extremely  high  degrees  of  interpolation  (p  —>  oo).  Substituting  the 
limiting  values  of  p  into  (11),  and  noting  that  equality  will  suffice  to  ensure  that  the 
required  accuracy  is  attained,  we  obtain  bounds  for  the  selection  of  N*,  namely 

A/.tt   <   N.   <   — 

for  all  values  of  p  >  0. 

The  existence  and  uniqueness  of  optimal  solutions  are  fairly  easy  to  establish. 
W(N*,p)  is  continuous  with  respect  to  each  of  its  variables,  and  both  of  the  first 
partial  derivatives  are  everywhere  positive.  This  observation  leads  to 

Lemma  2.  Let  S  be  the  set  {(N*,p)  :  X.  >  A/.7r(A701/(?>fl)}.  and  let  OS  be 
that  portion  of  the  boundary  of  S  given  by  <(X.,p)  :  X.  =  Mt,ir(N/e)1'^1^  >.    Then 

if(*o,yo)  G  S,  there  exists  a  point  (£.?/)  G  OS  such  that  JF(£.7/)  <  W(xo,yo). 

i 

Proof:  Since  (x0.y0)  G  S,  the  point  (£,y0)  G  OS.  where  £  =  Mmir  (Ey°+\ 
Furthermore,  f  <  xu.  Then  since  the  partial  derivative  of  the  work  function  with 
respect  to  Nm  is  everywhere  positive.  H"(£.yo)  <  U'(x0,yo)- 

The  utility  of  Lemma  2  is  that  the  optimization  problem  can  be  rewritten  as  a 
problem  in  a  single  variable.  Since  for  every  point  in  S  there  is  some  point  along  OS 
that  requires  less  work,  it  is  only  necessary  to  seek  a  minimum  from  the  points  of  OS . 
This  can  be  done  by  parameterizing  X.  and  p  as  functions  of  a  single  variable. 


Then  on  OS  we  find  that 
(14) 


.V.  =  M.nb       and       ;;+l=logfcf  —  J 


Since  0  <  p  <  oo,  the  value  of  b  is  restricted  to  the  interval  (1,  jV/e].  Substituting 
these  expressions  into  (12),  the  work  equation  may  be  rewritten  as  a  function  of  b 
alone 

(15)  \V(b)  =  Cu\L7ib\og2(Mm7ib)  +  2X\ogb(±)    , 

and  the  problem  is  to  minimize  (15)  subject  to  the  constraint  1  <  b  <  {X/e).  Once  b 
is  determined,  the  necessary  values  of  A',  and  p  can  be  obtained  from  (14).  We  may 
now  establish 

THEOREM   1.     There  exists  a  unique  value  b  that  minimizes  (15)  subject  to  1  < 
6  leq(X/e).   Therefore  the  work  function 

W(N„p)  =  d  Ar.  log2  X.  +  2X(p  +  1 ) 
has  a  unique  minimum,  subject  to  the  constraints 


Xm  >  A/.7T  I  —  J  Ofl(f     0<p<oo 
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Proof:   W{b)  is  continuous  and  difTerentiable  with  respect  to  b  on  (l,JV/c],  Dif- 
ferentiating equation  (15)  yields 

{U>)  \\"(b)  =  K1\n(K2b)- 


6(  In  6)2     ' 


where 


In  2  \  e 

For  W'(6)  =  0,  then.  6  must  satisfy  6(ln  6)2  ln(  A'26)  =  K$/ K\.  Now  W  is  also  contin- 
uous and  differentiate  on  (1, 7V/c],  and  differentiating  yields 


w«(h\      K]       r    flnb  +  2\ 


Since  6  >   1  we  see  that   \\'"(b)  >  0  for  all  b  £  (l,.V/(],  so  any  critical  point  in  the 
interval  must  correspond  to  a  local  minimum. 

It  is  apparent  that  U''(6)  — ►  -oc  as  b  —  1.  Examination  of  the  endpoint  b  —  N/e 
reveals  that  since  K\  >  ~,  I\2  >  e,  and  f<l.  we  have  that  \V'\  N/e)  >  0.  Further, 
\\'"(b)  >  0  implies  thai  W'{b)  has  exactly  one  sign  change  in  the  interval  {\,.\/(]. 
The  point  b0  at  which  this  occurs  is  therefore  a  global  minimum  for  Wib).  and  the 
value  lt'(  \..p).  where 

.V.  =  A/.7T&0      and      p  =  log,     (  —  )  -  1 


'»'>0     , 


is  the  unique  global  mini  mum  for  W  on  OS. 

The  values  of  .V.  and  p  obtained  in  this  manner  are  real  numbers.  There  is  only  a 
limited  number  of  integers  for  which  efficient  FFTs  exist,  and  Lagrangian  interpolation 
requires  p  to  be  an  integer.  Further,  this  entire  discussion  has  been  predicated  on  the 
assumption  that  p  is  an  odd  integer,  although  a  similar  analysis  can  be  made  for  p 
even.  Once  the  theoretical  values  of  N,  and  p  are  determined,  they  must  be  modified 
to  allow  computation.  There  is  some  flexibility  in  this,  but  certainly  selecting  N»  to 
be  the  first  integer  larger  than  M.irb  for  which  an  FFT  exists,  and  choosing  p  to  be 
the  smallest  odd  integer  greater  than 

N 

€ 


log,  (  -  -  )  -  1 


will  suffice. 

In  order  to  find  the  optimal  values  of  p  and  X.  it  is  necessary  to  find  the  value  of 
b  satisfying 

(IT)  b(\nb)2ln(K2b)  =  ^-  . 

While  an  analytic  solution  of  this  equation  cannot  be  found,  Newton's  iteration  may 
be  used.  Table  1  displays  optimal  parameters  Nm  and  p  for  several  combinations  of 
N,  A/,,  and  c. 


A 

A/. 

( 

A. 

P 

A' 

AL 

e 

A, 

P 

32 

8 

.1 

48.7 

7.7 

128 

32 

.1 

193 

9.9 

32 

32 

.1 

142.6 

15.5 

128 

64 

.1 

325.7 

13.8 

32 

64 

.1 

257.6 

22.2 

128 

128 

.1 

570.7 

19.4 

32 

8 

.01 

52.9 

9.8 

128 

32 

.01 

206.7 

12.1 

32 

32 

.01 

150.1 

19.1 

128 

64 

.01 

344.1 

16.6 

32 

64 

.01 

267.8 

27.1 

128 

128 

.01 

595.6 

23.1 

Table  1.  Optimal  parameters  Nm  and  p  computed  for  various  problems. 

5.  An  ADFT  Example.  To  illustrate  the  AD  FT,  consider  the  problem  of  com- 
puting the  formal  DFT  of  the  function  u(x)  —  [(w  —  x)/ir]2,  sampled  on  an  irregular 
grid.  The  irregular  grid  consists  of  A  =  128  points  x}  randomly  spaced  in  the  interval 
(0,27r).  Since  the  extent  of  the  interval  is  27T,  the  frequencies  u>/  are  just  the  integers 
/,  and  the  formal  DFT  is 


is: 


,V-1 

a(0  =  ]C  u(xj)e 

j  =  0 


-ilx, 


-64  <  /  <  64 


The  sampled  data  are  shown  in  the  top  of  Figure  1.  The  real  part  of  (18)  is  plotted 
on  the  bottom  of  Figure  1.  The  AD  FT  was  used  to  approximate  the  formal  DFT. 
with  values  of  A'.  =  128,  256.  512,  and  1024.  Figure  2  displays,  for  each  choice 
of  A'„,  the  absolute  value  of  the  error  \u(l)  —  [\Vu](l)\.  plotted  as  a  function  of  /. 
Linear  interpolation  (p  =  1)  was  used  in  each  case.  Note  that  increasing  the  value 
of  A',  produces  a  noticeable  decrease  in  the  error,  and  that  the  error  increases  with 
increasing  wavenumber,  as  might  be  inferred  from  (10).  Figure  3  displays  the  effect  of 
using  different  values  of  p  for  fixed  A„.  It  may  be  seen  that  the  error  decreases  rapidly 
as  p  is  increased.  Equation  (9)  predicts  that  the  error  should  decrease  at  least  as  fast 

as  ( -vr)  decreases  as  p  or  A',  are  increased.  Table  2  gives  both  the  infinity  and 
Li  norms  of  the  error  |u(/)  -  [U'u](/)|  for  several  values  of  p  and  each  of  A',  =  256 
and  N»  =  512.  For  Ar,  =  256.  the  error  bound  decreases  by  0.6169  each  time  p  is 
increased  by  2.  The  experimental  error  is  diminished  by  a  factor  of  approximately 
0.3  as  p  is  increased  from  1  to  3,  and  by  a  factor  of  approximately  0.4  with  each 
succeeding  increase,  better  than  the  theory  predicts.  Similarly,  for  N*  =  512,  the 
theoretical  bound  decreases  by  0.15421  as  p  is  increased  by  2.  while  the  experimental 
decrease  is  approximately  0.11  for  each  increase,  a  slightly  better  result.  Numerical 
experiments  on  numerous  other  irregularly  sampled  functions,  with  various  degrees 
of  smoothness,  produced  similar  results.  In  these  experiments  the  ADFT  behaved 
in  a  similar  fashion  as  it  did  for  the  function  discussed  above.  There  is  dramatic 
improvement  with  increasing  values  of  A'»,  and  p.  As  might  be  expected,  the  error 
diminished  faster  with  smooth  functions  than  discontinuous  functions. 


A. 


\Error\ 


\Error\\-2 


A. 


\Error\ 


Error]  |2 


256 
256 
256 
256 


1.20663 
0.357889 
0.144478 
0.061305 


0.434861 
0.116080 
0.039136 
0.014983 


512 
512 
512 
512 


0.290501 
0.029351 
0.003360 
0.000419 


0.109943 

0.008315 

0.000817 

9.66385e-05 


Table  2.  Errors  of  the  ADFT  for  various  values  of  N.  and  p. 
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6.  Some  Open  Questions  about  the  ADFT.  Like  the  continuous  Fourier 
transform,  the  DFT  has  several  important  properties,  such  as  linearity,  the  convolu- 
tion and  correlation  properties,  the  shifting  property,  the  modulation  property,  and 
Parseval's  relation.  To  what  extent  these  properties  hold  for  the  ADFT  is  an  open 
question.  The  linearity  holds  can  be  established  immediately,  by  noting  that  both  the 
formal  DFT  and  the  ADFT  can  be  written  as  matrix  operations,  so  they  are  linear 
operators.  Certain  symmetry  properties  are  easy  to  establish.  For  example,  applying 
the  ADFT  to  a  real-valued  vector  will  yield  a  conjugate  symmetric  result,  that  is 
ii(uj)  =  u(—  u>),  because  the  vector  [I*]  u  is  real-valued,  and  because  the  ADFT  is 
computed  by  applying  the  FFT  operator  to  this  vector.  The  DFT,  and  therefore  the 
FFT,  maps  a  real  vector  to  a  conjugate  symmetric  vector  [4].  Applying  the  DFT  to 
data  vectors  with  other  symmetries  (even,  odd,  quarter-wave,  etc.)  yields  output  vec- 
tors with  other  types  of  symmetries  [10].  It  is  natural  to  ask  which  of  these  symmetry 
properties  are  inherited  by  the  formal  DFT  or  the  ADFT.  It  seems  reasonable  to 
postulate  that  if  the  irregular  gridpoints  are  symmetrically  disposed  and  the  function 
u(Xj)  is  symmetric  then  the  symmetry  property  of  the  DFT  might  be  inherited  by 
thp  formal  DFT  and  the  ADFT. 

An  important  question  is:  How  is  the  formal  DFT  related  to  the  continuous 
Fourier  transform?  That  is,  to  what  extent .  and  with  what  prror.  does  the  formal  DFT 
approximate  the  FT1  Answering  this  question  may  prove  to  be  a  lengthy  process. 
Many  related  questions  will  also  arise  For  example,  how  does  the  sampling  theorem 
apply  to  an  irregular  grid?  What  frequencies  can  be  represented  accurately,  and  what 
constitutes  aliasing?  Is  there  some  analog  to  the  Poisson  summation  theorem?  Many 
problems  feature  irregularly  spaced  data,  so  it  may  be  assumed  that  these  questions 
are  of  some  interest. 
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THE  AXTZRPOLAXED  DFT 
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Fig.   1.   The  function  /(r)    =    [(x  -  x)/t]2  sampled  on  in  irregular  grid,  and  its  formal  DFT.  Only 
the  real  part  of  the  formal  DFT  is  plotted. 
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Fig.  2.    The  absolute  value  of  the  error  of  the  AD  FT.  plotted  as  a  function  of  wavenumber.    Each 
graph  corresponds  to  a  different  choice  of  .V.,  while  linear  interpolation  !>  =  1)  is  used  for  all. 
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Fig.  3.  The  absolute  value  of  the  error  of  the  ADFT,  plotted  as  a  function  of  wavenumber.  N,  =  512 
for  each  graph,  while  several  values  of  p  are  used. 
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